{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Medoid Compositing\n",
    "*Seasonal Composite Landsat TM/ETM+ Images Using the Medoid (a Multi-Dimensional Median), Neil Flood, 2013, doi:10.3390/rs5126481*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geetools import tools, composite, cloud_mask, indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipygee import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build a collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = ee.Geometry.Point(-72, -42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\\\n",
    "        .filterBounds(p).filterDate('2017-01-01', '2017-12-01')\\\n",
    "        .map(cloud_mask.landsat8SRPixelQA())\\\n",
    "        .map(lambda img: img.addBands(indices.ndvi(img,'B5', 'B4')))\\\n",
    "        .limit(7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ee460daf99e44237a1364eebcf16758a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Accordion(children=(Output(),), _titles={'0': 'Loading...'}),))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eprint(col.size())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other simple composites to compare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_ndvi = col.qualityMosaic('ndvi')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "mosaic = col.mosaic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Medoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### add date band before compositing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_date(img):\n",
    "    date = tools.date.getDateBand(img)\n",
    "    return img.addBands(date).copyProperties(date, ['day_since_epoch'])\n",
    "col = col.map(add_date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "minval = ee.Number(col.aggregate_min('day_since_epoch'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxval = ee.Number(col.aggregate_max('day_since_epoch'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "medoid = composite.medoid(col, bands=bands)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "medoid = medoid.set('min_date', minval, 'max_date', maxval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Medoid without taking in count zero values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "medoid_no_zeros = composite.medoid(col, bands=bands, discard_zeros=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "medoid_no_zeros = medoid_no_zeros.set('min_date', minval, 'max_date', maxval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show on Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1395e9b89dc148a5b3c7740862bf276b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(basemap={'max_zoom': 19, 'attribution': 'Map data (c) <a href=\"https://openstreetmap.org\">OpenStreetMap</a…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cefb9efff11f42d18fc299153049c35f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tab(children=(CustomInspector(children=(SelectMultiple(options=OrderedDict(), value=()), Accordion(selected_in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Map = Map()\n",
    "Map.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "vis = {'bands':['B5', 'B6','B4'], 'min':0, 'max':5000}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(p)\n",
    "Map.centerObject(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(max_ndvi, vis, 'max NDVI')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(mosaic, vis, 'simply Mosaic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(medoid, vis, 'Medoid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(medoid.select('date').randomVisualizer(), {'bands':['viz-red', 'viz-green', 'viz-blue'], 'min':0, 'max':255}, 'dates')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(medoid_no_zeros, vis, 'Medoid without zero values')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(medoid_no_zeros.select('date').randomVisualizer(), {'bands':['viz-red', 'viz-green', 'viz-blue'], 'min':0, 'max':255}, 'dates of medoid without zero values')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract data from images and compute locally to compare"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract medoid values in point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "medoid_values = tools.image.getValue(medoid.select(bands), p, scale=30, side='client')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'B2': 77, 'B3': 188, 'B4': 105, 'B5': 2500, 'B6': 691, 'B7': 224}"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "medoid_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List of values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "medoid_values_list = [val for _, val in medoid_values.items()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[224, 2500, 77, 691, 188, 105]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "medoid_values_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract values at point in each image of the collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_values = tools.imagecollection.getValues(col.select(bands), p, scale=30, side='client')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get bandnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_key_list = []\n",
    "for _, d in col_values.items():\n",
    "    keys = []\n",
    "    for k, v in d.items():\n",
    "        keys.append(k)        \n",
    "    col_key_list.append(keys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[['B7', 'B5', 'B2', 'B6', 'B3', 'B4'],\n",
       " ['B7', 'B5', 'B2', 'B6', 'B3', 'B4'],\n",
       " ['B7', 'B5', 'B2', 'B6', 'B3', 'B4'],\n",
       " ['B7', 'B5', 'B2', 'B6', 'B3', 'B4'],\n",
       " ['B7', 'B5', 'B2', 'B6', 'B3', 'B4'],\n",
       " ['B7', 'B5', 'B2', 'B6', 'B3', 'B4'],\n",
       " ['B7', 'B5', 'B2', 'B6', 'B3', 'B4']]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col_key_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get values as a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_values_list = []\n",
    "for _, d in col_values.items():\n",
    "    values = []\n",
    "    for _, v in d.items():\n",
    "        if v:\n",
    "            values.append(v)\n",
    "        else:\n",
    "            values.append(0)\n",
    "    col_values_list.append(values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[224.0, 2445.0, 83.0, 698.0, 200.0, 127.0],\n",
       " [287.0, 3168.0, 112.0, 870.0, 272.0, 143.0],\n",
       " [259.0, 2717.0, 90.0, 813.0, 210.0, 120.0],\n",
       " [224.0, 2500.0, 77.0, 691.0, 188.0, 105.0],\n",
       " [245.0, 2465.0, 87.0, 720.0, 193.0, 107.0],\n",
       " [0, 0, 0, 0, 0, 0],\n",
       " [307.0, 3142.0, 107.0, 928.0, 290.0, 159.0]]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col_values_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Medoid Method locally"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def local_medoid(values):\n",
    "    from copy import copy\n",
    "    import math\n",
    "\n",
    "    def distance(arr1, arr2):\n",
    "        zipped = zip(arr1, arr2)\n",
    "        accum = 0\n",
    "        for a, b in zipped:\n",
    "            calc = (a-b)*(a-b)\n",
    "            accum += calc\n",
    "        return math.sqrt(accum)\n",
    "\n",
    "    def med(values):\n",
    "        results = {}\n",
    "        for i, val in enumerate(values):\n",
    "            val = list(val)\n",
    "            cop = copy(values)\n",
    "            cop = [list(a) for a in cop]\n",
    "            cop.remove(val)\n",
    "            dist = 0\n",
    "            for r in cop:\n",
    "                r = list(r)\n",
    "                d = distance(val, r)\n",
    "                dist += d\n",
    "            results[i] = dist\n",
    "\n",
    "        return results\n",
    "    \n",
    "    def getmin(d):\n",
    "        minval = min(d.values())\n",
    "        for k, v in d.items():\n",
    "            if v == minval:\n",
    "                return k\n",
    "    \n",
    "    values = med(values)\n",
    "    min_value = getmin(values)\n",
    "    \n",
    "    # return the index of the minimized sum as first argument, and all options as second\n",
    "    return min_value, values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute medoid locally and compare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "local = local_medoid(col_values_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,\n",
       " {0: 4461.227704323096,\n",
       "  1: 6022.895734074118,\n",
       "  2: 4593.00466857401,\n",
       "  3: 4380.0310171598885,\n",
       "  4: 4399.422196724632,\n",
       "  5: 17251.127189606195,\n",
       "  6: 5996.407209928899})"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "local"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the values that correspond to the medoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_values = col_values_list[local[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[224.0, 2500.0, 77.0, 691.0, 188.0, 105.0]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Match bands with values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "local_medoid = dict(zip(col_key_list[0], min_values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'B2': 77.0, 'B3': 188.0, 'B4': 105.0, 'B5': 2500.0, 'B6': 691.0, 'B7': 224.0}"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "local_medoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'B2': 77, 'B3': 188, 'B4': 105, 'B5': 2500, 'B6': 691, 'B7': 224}"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "medoid_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finally, compare values from medoid mosaic against locally computed medoid (from images values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "medoid_values == local_medoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
